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Abstract 

We are concerned with high-fidelity subsurface imaging of the soil, which commonly 
arises in geotechnical site characterization and geophysical explorations. Specifically, we 
attempt to image the spatial distribution of the Lame parameters in semi-infinite, three- 
dimensional, arbitrarily heterogeneous formations, using surficial measurements of the soil’s 
response to probing elastic waves. We use the complete waveform response of the medium 
to derive the inverse problem, by using a partial-differential-equation (PDE)-constrained 
optimization approach, directly in the time-domain, to minimize the misfit between the 
observed response of the medium at select measurement locations, and a computed response 
corresponding to a trial distribution of the Lame parameters. We discuss strategies that 
lend algorithmic robustness to our proposed inversion scheme. To limit the computational 
domain to the size of interest, we employ perfectly-matched-layers (PMLs). 

In order to resolve the forward problem, we use a recently developed hybrid finite 
element approach, where a displacement-stress formulation for the PML is coupled to a 
standard displacement-only formulation for the interior domain, thus leading to a com¬ 
putationally cost-efficient scheme. Time-integration is accomplished by using an explicit 
Runge-Kutta scheme, which is well-suited for large-scale problems on parallel computers. 

We verify the accuracy of the material gradients obtained via our proposed scheme, 
and report numerical results demonstrating successful reconstruction of the two Lame 
parameters for both smooth and sharp profiles. 
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1. Introduction 


Seismic inversion refers to the process of identification of material properties in geo¬ 
logical formations [8, 31, 36]. The problem arises predominantly in exploration geophysics 
[14, 23, 25, 33] and geotechnical site characterization [20]; it belongs to the broader class of 
inverse medium problems: waves, whether of acoustic, elastic, or electromagnetic nature, 
are used to interrogate a medium, and the medium’s response to the probing is subse¬ 
quently used to image the spatial distribution of properties (e.g.. Lame parameters, or 
wave velocities) [4, 15, 21]. Mathematically, algorithmically, and computationally, inverse 
medium problems are challenging, especially, when no a priori constraining assumption 
is made on the spatial variability of the medium’s properties. The challenges are further 
compounded when the underlying physics is time-dependent, and involves more than a 
single distributed parameter to be inverted for, as in seismic inversion. 

Due to the complexity of the inverse problem at hand, most techniques to date rely 
on simplifying assumptions, aiming at rendering a solution to the problem more tractable. 
These assumptions can be divided into five categories: a) assumptions regarding the di¬ 
mensionality of the problem, whereby the original problem is reduced to a two-dimensional 
[17, 20, 22], or a one-dimensional problem [26]; b) assuming that the dominant portion of 
the wave energy on the ground surface is transported through Rayleigh waves, and thus, 
disregarding other wave types, such as compressional and shear waves, as is the case in the 
Spectral-Analysis-of-Surface-Waves (SASW) and its variants (MASW) [35]; c) inverting 
for only one parameter, as is done in [1, 11, 28, 29], where inversion was attempted only 
for the shear wave velocity, assuming the compressional wave velocity (or an equivalent 
counterpart) is known; d) assumptions concerning the truncation boundaries, which are 
oftentimes, grossly simplified due to the complexity associated with the rigorous treatment 
of these boundaries [40]; and e) idealizing the soil body, which is a porous and lossy medium, 
as an elastic solid and neglecting its attenuative properties^ [12]. Over the past decade, 
continued advances in both algorithms and computer architectures have allowed the grad¬ 
ual removal of the limitations of existing methodologies. However, a robust methodology, 
especially for the time-dependent elastic case remains, by and large, elusive. 

Among the recent works on inversion, which are similar in character to ours, we refer 
to Pratt et al. [32] who considered two-dimensional acoustic inversion in the frequency 
domain, and Epanomeritakis et al. [11] where full-waveform inversion has attempted for 
three-dimensional time-domain elastodynamics, where a simple boundary condition was 
used for domain truncation. Kang and Kallivokas [21] considered the problem for the two- 
dimensional time-domain acoustic case, and used PMLs to accurately account for domain 
truncation. Kucukcoban [22] extended the work of Kang and Kallivokas to two-dimensional 
elastodynamics, and reported successful reconstruction of the two Lame parameters for 


^See [2] for a full-waveform-inversion-based approach, using a generalized Maxwell model for lossy soils 
in a one-dimensional setting. 
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models involving synthetic data. Bramwell [7] used a discontinuous Petrov-Galerkin (DPG) 
method in the frequency domain, endowed with PMLs, for seismic tomography problems, 
advocating the DPG scheme over conventional continuous Galerkin methods, since it re¬ 
sults in less numerical pollution. Recently, Jung et al. [18, 19] used an extended finite 
element method (XFEM) to explicitly parameterize the boundaries of scatterers for two- 
dimensional problems in elastodynamics. Their approach seems promising especially for 
the identification of voids. 

In this paper, we discuss a systematic framework for the numerical resolution of the 
inverse medium problem, directly in the time-domain, in the context of geotechnical site 
characterization. The goal is to image the arbitrarily heterogeneous material profile of 
a probed soil medium, using complete waveforms^ of its response to interrogating elastic 
waves, originating from the ground surface. To this end, the response of the soil medium 
to active sources (Vibroseis equipment) is collected by receivers (geophones) dispersed over 
the formation’s surface, as shown in Figure 1(a). Arriving at a material profile is then 
accomplished by minimizing the difference between the collected response at receiver lo¬ 
cations, and a computed response corresponding to a trial distribution of the material 
parameters. Imaging near-surface deposits brings additional difficulties, typically not en¬ 
countered in exploration geophysics. In geophysical explorations, the probing is over large 
length scales; thus, an accurate domain termination tool may not play a critical role. 
However, in geotechnical site characterization, one, typically, deals with a much smaller 
domain. Moreover, obtaining a high-fidelity image of the near-surface deposits has practi¬ 
cal significance for the safe founding of infrastructure components; thus, having accurate 
termination conditions becomes critical. In this vein, and in the presence of heterogeneity, 
using PMLs for domain termination is the best available option, and is thus used in this 
work. The PML is a buffer zone that surrounds the domain of interest, and enforces the 
decay of outgoing waves. Figure 1(b) shows a computational model, obtained through the 
introduction of PMLs at the truncation boundaries. 

In order to address all the difficulties outlined earlier, we integrate recent advances in 
several areas. Specifically, we use (a) a recently developed state-of-the-art parallel wave sim¬ 
ulation tool for domains terminated by PMLs, which renders the computational model asso¬ 
ciated with the near surface geotechnical investigations finite [13]; (b) a partial-differential- 
equation (PDE)-constrained optimization framework through which the minimization of 
the difference between the collected response at receiver locations and a computed response 
corresponding to a trial distribution of the material properties is achieved [30]; (c) regular¬ 
ization schemes to alleviate the ill-posedness inherent in inverse problems; (d) continuation 
schemes that lend algorithmic robustness [21]; and (e) a biasing scheme that accelerates 
the convergence of the A-profile for robust simultaneous inversion of both Lame parameters 


^Using the complete waveform (complete recorded response) results in a full-waveform inversion ap¬ 
proach. 
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Figure 1: Problem definition: (a) interrogation of a heterogeneous semi-infinite domain by an active source; 
and (b) computational model truncated from the semi-infinite medium via the introduction of PMLs. 


[ 22 ]. 

The remainder of this paper is organized as follows: first, we discuss the numerical 
resolution of the forward problem. Next, we discuss the mathematical and numerical 
aspects of the underlying inverse medium problem, where we derive the adjoint and control 
problems, and discuss strategies that invite robustness. We then verify the accuracy of the 
material gradients that we compute, by comparing them with directional hnite differences. 
We report on numerical experiments, using synthetic data, targeting the reconstruction of 
both smooth and sharp profiles. Lastly, we conclude with summary remarks. 

2. Forward wave simulation in a 3D PML-truncated medium 

In the forward problem, we are concerned with the propagation of elastic waves in a 
three-dimensional, semi-infinite, arbitrarily heterogeneous elastic medium. To keep the 
computation feasible, one needs to limit the extent of the computational domain. This 
can be accomplished by placing PMLs at the truncation boundaries. Theoretically, the 
truncation boundary is reflectionless, and when outgoing waves pass through the interface, 
they get attenuated within the PML buffer zone. This concept is schematically captured 
in Figure 2. 

For the numerical resolution of the forward problem, we use a recently developed hybrid 
approach that uses a displacement-stress formulation for the PML buffer, coupled with a 
standard displacement-only formulation in the interior domain^. This approach results in 


3 


The terms “interior domain” and “regular domain” are used interchangeably throughout this article. 
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Figure 2: A PML truncation boundary in the direction of coordinate s. 


a computationally cost-efficient scheme, due to the treatment of the interior domain with 
a standard displacement-only elastodynamics formulation. We refer to [13] and references 
therein for the complete development and parallel implementation of the method. Herein, 
we only outline the resulting coupled system of equations. 



Figure 3: PML-truncated semi-infinite domain. 


Accordingly, find u(x, t) in U and S(x, t) in (see Figure 3 for domain 

and boundary designations), where u and S reside in appropriate function spaces and: 


5 















div {/i [Vu + (Vu)^] + A(divu)X} + b = pii in x J, 

(la) 

div ^S^Ae + S^Ap + S^A^^ = p (aii + 5u + cu + du) in x J, 

(lb) 

fflS + 5S + cS + dS = 

p [(Vu)Ae + Ae(Vu)^ + (Vu)Ap + Ap(Vu)^ + (Vu)A^ + Au,(Vu)^] + 

A [div(Aeu) + div(Apu) + div(Au,u)] I in 11™^ x J. 

(l c) 

The system is initially at rest, and subject to the following boundary and interface condi¬ 
tions: 


{/i [Vu -h (Vu)^] -|- A(div u)X} n 

9n 

on P^^ X J, 

(2a) 

(S^Ae -f S^Ap -f S^A^)n = 0 


on P^“L X J, 

(2b) 

u = 0 


on Pg“L X J, 

(2c) 

^RD _ ^PML 


on P^ X J, 

(2d) 

{/r [Vu -1- (Vu)^] + A(div u)X} n 

= (S^Ae + S^Ap + S^A^)n 

on P^ X J, 

(2e) 


where temporal and spatial dependencies are suppressed for brevity; u is the displacement 
vector, p denotes mass density of the medium, A and p are the two Lame parameters, 
X is the second-order identity tensor, S represents the Cauchy stress tensor, a dot (') 
denotes differentiation with respect to time, and a bar (“) indicates history of the subtended 
variable"^; denotes the interior (regular) domain, represents the region occupied 

by the PML buffer zone, L^ is the interface boundary between the interior and PML 
domains, P^^ and P™^ denote the free (top surface) boundary of the interior domain 
and PML, respectively, J = (0,r] is the time interval of interest, b denotes body force per 
unit volume, and is the prescribed surface traction. Moreover, Ae, Ap, and A^^, are the 
so-called stretch tensors, which enforce dissipation of waves in and a, b, c, and d are 

products of certain elements of the stretch tensors [13]. Eq. (la) is the governing PDE for 
the interior elastodynamic problem, whereas Eqs. (lb) and (Ic) are the equilibrium and 
combined kinematic and constitutive equations, respectively, for the PML buffer. 

Next, we seek a weak solution, corresponding to the strong form of (1) and (2), in the 
Galerkin sense. Specifically, we take the inner products of (la) and (lb) with (vector) test 


"^For instance, u(x,t) = u(x, t) dr. 
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function w(x), and integrate by parts over their corresponding domains. Incorporating 
(2d-2e) eliminates the interface boundary terms and results in (3a). Next, we take the 
inner product of (Ic) with (tensor) test function T(x); there results (3b). We refer to [13] 
for details. 

Accordingly, find u G H^(fl) x J, and S G £^(11) x J, such that: 


J Vw : {/r [Vu + (Vu)^] + A(divu)X} dfl + J Vw : ^S^Ae + S^Ap + S^A^,^ dfl 

QRD QPML 

+ / w ■ pii dfl + w ■ p (aii + 6u + cu + du) dfl = w ■ dT + /w-b dfl, (3a) 


nRD 


QPML 




nRD 


j f-.(^aS + bS + cS + dS^ dQ 

QPML 

= j T :p [(Vu)Ae + Ae(Vu)^ + (Vu)Ap + Ap(Vu)^ + (Vu)A^ + A^(Vu)'^] 

qPML 

+ T : A [div(Aeu) + div(Apu) + div(A^u)] I dfl, (3b) 

for every w G H^(n) and T G £^(n), where G L^(n) x J, and b G L^(n) x J. Function 
spaces for scalar-valued (u), vector-valued (v), and tensor-valued (^) functions are defined 
as: 


L^iyi) = |u: j jupdx < oo| , (4a) 

l2(II) = {v: vG(x2(fI))3}, (4b) 

£2(0) = {^: ^ G (£2(0))3x3} , (4c) 

Fr^(fl) = |u: j (|u|2-|-|Vt>|2) dx < oo, u|pPML= o| , (4d) 

= {v: V e . (4e) 


In order to resolve (3) numerically, we use standard finite-dimensional subspaces. Specif¬ 
ically, we introduce finite-dimensional subspaces c H^(0) and c £2(0), with basis 
$ = (<I>i, <I' 2 ,..., 4>Ar)^ and ^ = ('hi, '1^2,..., 'I'm)^, respectively, where N and M denote 
the dimension of of the corresponding vector space. We then approximate u(x, t) with 
u/j(x, t) G X J, and S(x, t) with S/i(x, t) G ~fh x J. In a similar fashion, we approximate 
the test functions, w(x) with w/j(x) G and T(x) with T/i(x) G T/^. It can be shown 
[13] that the following semi-discrete form results: 
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(5a) 

(5b) 


Md"* + Cd"* + Kd"* + Gd"* = P\ 

d®* = f d®*(T)|pML d'T, 

Jo 

where spatial and temporal dependencies are suppressed for brevity; M, C, K, G, are 
system matrices, d®* = (u^, S^)^ is the vector of nodal unknowns comprising displacements 
in and stress components only in and P* is the vector of applied forces. 

The matrix M has a block-diagonal structure, and can be diagonalized if one employs 
spectral elements with a Legendre-Gauss-Lobatto (LGL) quadrature rule, which then en¬ 
ables explicit time integration of (5). In this regard, we express (5) as a first-order system: 


Xo ' 


■ 0 

I 

0 ■ 


Xo' 


o' 

Xl 

= 

0 

0 

I 


Xl 

+ 

0 

Mx2 


-G 

-K 

-c 


X2 


^st 


where xq = d®*, xi = d®*, and X 2 = d®*. We then use an explicit fourth-order Runge-Kutta 
(RK-4) method for integrating (6) in time. 

3. The inverse medium problem 

Our goal is to find the distribution of the Lame parameters A(x) and fi(x) within 
the elastic soil medium. We consider sources, and the response recorded at receivers on 
the ground surface, as known. The inverse medium problem can thus be formulated as 
the minimization of the difference (or misfit) between the measured response at receiver 
locations, and a computed response corresponding to a trial distribution of the material 
parameters. The misfit minimization should honor the physics of the problem, which is 
idealized by the forward problem, stated in the preceding section. Mathematically, this 
can be cast as a PDE-constrained optimization problem: 

1 Wr „Tr- 

min J(A,/r) := 77 V] / / {u - Um) ■ {n - Um) S{x - x^) dT dt + n{X, fi), (7) 

where u is the solution of the forward problem governed by the initial- and boundary-value 
problem (1), (2). 

In the above, J is the objective functional^, W denotes the total number of receivers, 
T is the total observation time, T^ is the part of the ground surface where the receiver 


®We use J to indicate the corresponding discrete objective functional. See [6, 10] for other possibilities. 











response, u^, has been recorded, 5{x — x^-) is the Dirac delta function, which enables 
measurements at receiver locations Xj, and TZ{X,^) is the regularization term, which is 
discussed below. 

Inverse problems suffer from solution multiplicity, which, in general, is due to the 
presence of insufficient data. This makes the problem ill-posed in the Hadamard sense. 
Regularization of the solution by using the Tikhonov (TN) [37], or, the Total Variation 
(TV) [34] scheme are among common strategies to alleviate ill-posedness. The Tikhonov 
regularization, denoted by {X, fj,), penalizes large material gradients and, thus, pre¬ 
cludes spatially rapid material variations from becoming solutions to the inverse medium 
problem. It is defined as: 

jvxvxdn + 

where Rx and are the so-called A- and /i-regularization factor, respectively, and control 
the amount of penalty imposed via (8) on the gradients of A and /r. By construction, TN 
regularization results in material reconstructions with smooth variations. Consequently, 
sharp interfaces may not be captured well when using the TN scheme. The TV regulariza¬ 
tion, however, works better for imaging prohles involving sharp interfaces, as it typically 
preserves edges. It is defined as: 

7^^^(A,M) = ^ y’(VA-VA + e)5dD + ^ j{Vfx ■ V+ e)^ dQ, (9) 

ORD QRD 

where the parameter e makes differentiable when either VA ■ VA = 0, or, V/r ■ V/r = 0. 

For computing the first-order optimality conditions for (7), we use the (formal) La- 
grangian approach [41] to impose the PDE-constraint in its weak form. These are necessary 
conditions that must be satisfied at a solution of (7). Specifically, we introduce Lagrange 
multiplier vector function w G H^(D), and Lagrange multiplier tensor function T G 
to enforce the initial- and boundary-value problem (1), (2), which admits the weak form 
(3). The Lagrangian functional becomes: 


K™(A,rt = :| 


Rn 


■ V/r dD, 


( 8 ) 
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1 Wr „T ^ 

£(u,S,w,T, A,m) = 2 X] / / (u - Um) ■ (u - u^) 5(x - Xj) dr dt + 7^(A,//) 

j = l ^ 


rT 


rT 


j Vw; {/It [Vu + (Vu)^] + A(divu)Z} dfJ dt 

RD 

J Vw: 


S^Ae + S^Ap + S^A^ ) dndt- 


w ■ pii dO dt 


QPML 




J J w ■ p (aii + &u + cu + du) dfi dt + y J w ■ g^dT 


dt 


QPML 


pRD 
^ N 


rT 


+ / w ■ b dfi dt - / / T : aS + 6S + cS + dS dn dt 


QR-D 


/O 


QPML 


+ 


y J T -.p [(VA)A, + A,(Vii)^ + (Vu)Ap + Ap(Vu)^ + (Vu)A^ + A^(Vu)^] 

(10a) 

where now u, S, A, and p are treated as independent variables. 


QPML 

+ T : A [div(Aeu) + div(Apu) + div(Au,u)] I dfl dt, 


3.1. Optimality system 

We now use the Lagrangian functional (10) as a tool to compute the optimality system 
for (7). To this end, the Gateaux derivative® (or first variation) of the Lagrangian functional 
with respect to all variables must vanish. This process is discussed next. 


3.1.1. The state problem 

Taking the derivatives of the Lagrangian functional C with respect to w and T in 
directions w G H^(n) and T G and setting it to zero, results in the state problem, 

which is identical to (3). That is: 

£'(u,S,w,T,A,/r)(w,T) = 0. (11) 

3.1.2. The adjoint problem 

We now take the derivative of C with respect to u and S in directions u G H^(n) and 
S G £2(0). This yields: 


®See Appendix A for the definition and notation. 
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u-(u — Um) — Xj) dr dt 


Nr I.X 

£'(u,S,w,T,A,/i)(u,S) = ^ 

i=i 

rT 


0 JVr, 


j Vw : [Vu + (Vu)^] + A(div u)l} dfJ dt 

f Vw: (S 

^ML 

J J w ■ p ^au + 5u + cu + dvL^ 


QPML 


S^Ae + S^Ap + S^A^ ) dr! dt - 


w ■ pu df! dt 


QRD 


dr! dt 


QPML 


r 


aS + bS + cS + dS] dr! dt 


QPML 


rT 


+ 


T:m 


(Vu)Ae + Ae(Vu)^ + (Vu)Ap + Ap(Vu)-' + (Vu)Au; + A^(Vu) 


~\T 


QPML 


T:A 


div(Aeu) + div(Apu) + div(Au;u) 


X dr! dt. 


(12a) 


Setting the above derivative to zero, and performing integration by parts in time, results 
in the statement of the weak form of the adjoint problem. That is, find w G H^(r!) x J, 
and T e £^(r!) x J, such that: 
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Vu: {/i [Vw + (Vw)^] + A(div w)X} dn+ 




+ y u ■ pw dn + J u - p (aw — bw + cw — dw) dil 

qRD qPML 

~iT \ I Ti A I rr-iT A rp A '^T 


Vu: p 


-TAe - Ae + TAp + Ap - TA^ - A^ 


qPML 

+ A 

N, 


— T : div(Aeu) + T : div(Apu) — T : div(A^u) 

= / u-(u - Um) 5(x - Xj) dr, 

=1 


idn 


(13a) 


j Vw:S^Ae 

QPML 


Vw : S^Ap + Vw; S^A^ dQ = j S : (^aT - bt + cS - dx) dQ 

QPMh 

(13b) 


for every u G H^(r2) and S G £^(fl), where w(x, T) = 0, and T(x, T) = 0. 

We remark that the adjoint problem (13) is a final-value problem and, thus, is solved 
backwards in time^; it is driven by the misfit between a computed response, and the 
measured response at receiver locations. Moreover, the operators implicated in the adjoint 
problem are very similar to those of the state problem: they involve transposition of 
the system matrices, and sign reversal for terms involving history, and first-order time 
derivatives. In this regard, we obtain the following semi-discrete form for the adjoint 
problem: 


Mciadj _ (14a) 

rdadj(r)|pMLdT, (14b) 

Jo 

where superscript “adj” refers to the adjoint problem, = (w^,T^)^ is the vector of 
nodal unknowns comprising discrete values of w in 14^^ U and discrete values of T 

only in and is a vector comprising the misfit at receiver’s locations. Moreover, 

system matrices M, C, K, G, are identical to those of the forward problem and, thus. 


^See [44] for other possibilities, and [24, 38, 39] for alternative approaches. 
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with minor adjustments, an implementation of the forward problem can also be used for 
the solution of the adjoint problem. 

The matrix M in (14) can be diagonalized by using spectral elements with a Legendre- 
Gauss-Lobatto (LGL) quadrature rule, similar to what we did in (6). We rewrite (14) as 
a first-order system: 


d 

dt 


’ yo ' 


yi 

= 

_My2_ 

_ 

d^dj. 

y2 = 


0 10 

0 0 1 

qT _^t 



Vo' 


■ 0 ■ 


yi 

y2. 

-f 

0 

fadj 


(15) 


where yo = yi = y 2 = d^“b with final values yo(T') = 0, yi(T) = 0, and 

y2(r) = 0. We use an explicit RK-4 method to integrate (15) in time. The scheme is 
outlined in Appendix B. 


3.1.3. The control problems 

Lastly, we take the derivative of £ with respect to A and fj, in directions A and fl, which 
yields the reduced gradients with respect to A and fj,, respectively. We restrict the reduced 
gradients to (The material properties at the interfaces T^ are extended into the PML 
buffer.). For the TN regularization, this yields: 


£^(u, S, w, T, A,/i)(A) = i?A yVA-VAdri —y y A Vw :(div u)X dfl dt, (16a) 

nan ^ QRD 

cT 

a . 

/o 

QRD QRD 

(16b) 


2'(u, S, w, T, A, /i)(/i) = J Vjl ■ Vp dn — J j jl Vu : [Vw -1- (Vw)^] dll dt. 


Setting the above derivatives to zero, results in the control problems. Similarly, for the TV 
regularization, the control problems read: 


£'(u,S,w,T,A,m)(A) = i?A 

£'(u, S, w, T, A, = Rf, 


VA ■ VA f ~ 

-^ dn — / / A Vw :(div u)Z dfl dt, (17a) 

J (VA-VA-Fe )2 do J 

ORD ^ ^ qRD 


V/i ■ V/x 


rT 




(V/x-V/x + e )2 Jo 


- dn — jl Vu: [Vw + (Vw) ] dn dt. 




(17b) 


Discretization of either (16) or (17) result in the following form: 
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Mg^ = Rx greg + gmis> ( 18 a) 

Mg'^ = Rf, g(4g + (18b) 

where M is a mass-like matrix, g^ and g^ is the vector of discrete values of the (reduced) 
gradient for A and /u, respectively, and g^gg, greg and g^jg, g((,is are the associated vectors 
corresponding to the regularization-part and misfit-part of g^ and g'^. We refer to Appendix 
C for matrix and vector definitions, and discretization details. 

3.2. The inversion process 

A solution of (7) requires simultaneous satisfaction of the state problem (6), the adjoint 
problem (15), and the control problems (18). This approach -a full-space method- is, in 
principle, possible [9]; however, the associated computational cost can be substantial. Al¬ 
ternatively, a reduced-space method may be adopted, in which, discrete material properties 
are updated iteratively, using a gradient-based minimization scheme. The latter approach 
is employed here, and is discussed next. 

We start with an assumed initial spatial distribution of the control parameters (A and 
fj,), and solve the state problem (6) to obtain d®* = (u^, S^)^. With the misfit known, we 
solve the adjoint problem (15) and obtain = (w^,T^)^. With and known, the 
(reduced) material gradients, i.e., g^ and g^, can be computed from (18). Thus, the vector 
of material values, at iteration k + 1, can be computed by using a search direction via: 


Afc+i — Afc -|- s^, (19a) 

fJ-k+i = Mfc + (19b) 

where A and comprise the vector of discrete values for A and fj,, respectively, a^, are 
step lengths, and s^, are the search directions for A^ and Herein, we use the L-BFGS 
method to compute the search directions [27]®. Moreover, to ensure sufficient decrease of 
the objective functional at each inversion iteration, we employ an Armijo backtracking line 
search [27], which is outlined in Algorithm 1. The inversion process discussed thus far is 
summarized in Algorithm 2. 

We remark that for the reduced-space method, either (16) or (17) can also be expressed 
as: 


£'(u, S, w, T, A, = J'(A, /r)(A), (20a) 

£'(u, S, w, T, A, = J'{\ (20b) 


®In the numerical experiments that we perform, we store m = 15 L-BFGS vectors. 
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Algorithm 1 Backtracking line search. 

1: Choose a^, ci, p > e.g., = 1, = 1, ci = 10“'^, p = 0.5 

2: while J(Afc + fJ-k + Mfc) + ci(a'^ gfc ' Sfc ' 

3: pa^ 

4: <— pa^ 

5: end while 

6: Terminate with 


Algorithm 2 Inversion for Lame parameters. 

1 

A: ^ 0 


2 

Set initial guess for material property vectors A^, /i^ 


3 

Compute J(Afc,/r^) 

>Eq. (7) 

4 

Set convergence tolerance tol 


5 

while {J(Afc,/i^) > tol} do 


6 

Solve the state problem for d®* = {uj}, 

> Eq. (6) 

7 

Solve the adjoint problem for d®''^-' = (w^,T^)^ 

i> Eq. (15) 

8 

Evaluate the discrete reduced gradients g^, g^ 

> Eqs. (18) 

9 

Compute search directions s^, s}} 

0 L-BFGS 

10 

Choose step lengths a^, 

> Algorithm 1 

11 

Update material property vectors A^, 

i> Eq. (19) 

12 

A: A: + 1 


13 

end while 
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where the equality in (20) is due to the satisfaction of the state problem. Therefore, the 
reduced gradients in (18), are, indeed, the gradients of the objective functional with respect 
to A and fj,. 

3.3. Buttressing schemes 

Inverse medium problems are notoriously ill-posed. They suffer from solution multi¬ 
plicity; that is, material profiles that are very different from each other, and, potentially 
non-physical, can become solutions to the misfit minimization problem. Regularization 
of the control parameters alleviates the ill-posedness; however, this alone, may not be 
adequate when dealing with large-scale complex problems. In this part, we discuss addi¬ 
tional strategies that further assist the inversion process, outlined in Algorithm 2, to image 
complex profiles. 

3.3.1. Regularization factor selection and continuation 

Computation of the (reduced) gradients (18) necessitates selection of the regularization 
factors R\ and R^. A common strategy is to use Morozov’s discrepancy principle [42], 
where a constant value for the regularization factor is used throughout the inversion process. 
Here, we discuss a simple and practical approach that was initially developed for acoustic 
inversion [21], and, later, was successfully applied to problems involving elastic inversion 
[ 22 ]. 

We start by rewriting the discrete control problem (18), either for A or in the following 
generic form: 


Mg = R greg gmis, (21) 

where g refers to the vector of discrete values of the (reduced) gradient, either for A or 
M, greg and gmis are the associated vectors corresponding to the regularization-part and 
misfit-part of g, and R is the regularization factor yet to be determined. The main idea is 
that the “size” of R greg should be proportional to that of gmis at each inversion iteration. 
We define the following unit vectors for the two components of the gradient vector: 


where 


^^reg — 


greg 

llgregll’ 


^ _ Smis 

^mis — Ti iT ’ 
II Smis II 


denotes the Euclidean norm. Equation (21) can then be written as: 


( 22 ) 


Mg 


R IlSregll I^reg “ 1 “ ||Smis|| ^mis 

(23a) 

II Smis II II II ^reg ^misl 


^ II Smis II 


IlSmisll fp ^^reg “1“ ^misj 5 

(23b) 
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where, 


p = R (23c) 

II Smis II 

In (23b), for the “size” of p rij-eg to be proportional to ti mis throughout the entire inversion 
process, one may choose 0 < p < 1. Once a value for p has been decided on, R can be 
computed via: 


= (24) 

IlSregll 

where p can take large values (e.g., 0.5)® at early stages of inversion and, thus, narrow 
down the initial search space. As inversion evolves, p can be continuously reduced (e.g., 
down to 0.3) to allow for profile refinement. The suggested values for p are based on our 
experience with various numerical experiments that provide quality solutions for different 
test cases, and seem to be independent of the dimensionality, size, and discretization of the 
considered cases. 

3.3.2. Source-frequency continuation 

Using loading sources with low-frequency content result in an overall image of the 
medium that lacks fine features. To allow for more details, and fine-tune the profile, one 
needs to use sources with higher frequency content. Thus, the inversion process can be 
initiated with a signal having a low-frequency content and, then, the frequency range can be 
increased progressively as inversion evolves. This can be achieved by using a set of probing 
signals, ordered such that each signal has a broader range of frequencies than the previous 
ones. The inversion process then begins with using the first signal. Upon convergence, the 
prohle is used as a starting point with the second signal, and the process is repeated for 
all signals. 

3.3.3. Biased search direction for A 

Simultaneous inversion for both A and /r is remarkably challenging [11]. As we demon¬ 
strate in Section 4.2, the objective functional (7) is more sensitive to /r, than to A. Con¬ 
sequently, as the inversion evolves, the /r-profile converges faster than that of A. In [22], a 
biasing scheme was proposed to accelerate the convergence of the A-profile, such that, at 
the early stages of inversion, the search direction for A is biased according to that of /r. 

The main idea is that due to physical considerations, the A-profile should be, more 
or less, similar to the p-profile. Hence, during the early inversion iterations, the search 
direction for A is biased according to: 


®Since we normalize the regularization-part and misfit-part of the gradient in (23b), p = 0.5 means the 
gradient is weighted twice by the misfit than the regularization. 
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( 25 ) 


^ iPfcll iPfclr 

where W is a, weight that imposes the biasing amount. We assign full weight (W = 1) 
on fj, at the first inversion iteration, and reduce it linearly down to zero as iterates evolve 
(say at k = 50). After that, we let A evolve on its own, according to the original, unbiased 
search direction. 

4. Numerical experiments 

We present numerical experiments^^ with increasing complexity to test the proposed 
inversion scheme. In the first example, we verify the accuracy of the gradients, computed 
by using Algorithm 2. Next, we focus on material profile reconstruction for heterogeneous 
hosts, using synthetic data at measurement locations. Specifically, we consider: (a) a 
medium with smoothly varying material properties along depth, to study various aspects 
of the inversion scheme; (b) a horizontally layered profile with sharp interfaces; (c) a 
horizontally layered profile with an ellipsoidal inclusion, using highly noisy data; and (d) 
a layered profile with three inclusions in an attempt to implicate arbitrary heterogeneity. 
Throughout, we use Gaussian pulses to probe the considered domains: 


where the parameters that characterize the load are given in Table 1; /I is the mean, a 
is the deviation, fmax is the maximal frequency content of the pulse, tend is the active 
duration of the Gaussian pulse, and the load has an amplitude of 1 kPa. The time history 
of the loads and their corresponding Fourier spectrum are shown in Figure 4. 


Table 1: Characterization of Gaussian pulses. 


Name 

fmax 

P 

(7 

^end 

P2Q 

20 

0.11 

0.0014 

0.20 

P30 

30 

0.08 

0.0007 

0.15 

Pw 

40 

0.06 

0.0004 

0.12 


4 . 1 . Numerical verification of the material gradients 

Accurate computation of the discrete gradients is crucial for the robustness of Algo¬ 
rithm 2. The gradients of the objective functional with respect to the control parameters 
can be computed either by the optimize-then-discretize, or, the discretize-then-optimize ap¬ 
proach [30]. While the discretize-then-optimize method yields the exact discrete gradients 


developed a code in Fortran, using PETSc [3] to facilitate parallel implementation. 
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Figure 4: Time history of the Gaussian pulses and their Fourier spectrum. 


of the discrete objective functional [16], this is not always the case with the optimize-then- 
discretize scheme [43]. 

In this part, through a numerical experiment, we demonstrate that the discrete gradi¬ 
ents that we compute via the optimize-then-discretize technique, are accurate, and equal 
to the discrete gradients of the discrete objective functional. To this end, we compare 
directional hnite differences of the discrete objective functional, with directional gradients 
obtained from (18). We start by defining the hnite difference directional derivatives: 


where A and p, is the discrete direction vector for A and /i, respectively. The directional 
derivatives obtained via the control problems (18) are: 

d™(A, = A^ M g^, (27a) 

= Mgr (27b) 


Next, we verify that (26) and (27) produce identical values for several choices that we 
make for A and /i, by considering a test problem displayed in Figure 5, and detailed below, 
with regularization factors Rx = = 0^^. We consider perturbations A or p: the unit 


11 


Zero values are considered since convergence difficulties that may arise stem from the misfit part of the 
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vector is zero everywhere except at the component corresponding to coordinate {x, y, z) 
where the directional derivatives are being computed. The derivatives and with 
respect to either A or y,, for points with coordinates {x,y,z), are presented in Table 2. 
Digits where df^ coincides with are shown in bold. Since pointwise perturbations result 
in small changes in the objective functional, numerical roundoff influences the accuracy of 
the finite difference directional derivatives, as it has also been reported in [43]. 



Figure 5: Problem configuration for the verification of the gradients. 


Table 2: Comparison of the directional derivatives. 


Case 

fmax 

(x,y,z) 

Pert. 

field 


1 

o 

II 

h = 10"'‘ 

h = 10“® 

1 

20 Hz 

(1,1.0) 

A 

-3.03500e-9 

-3.03416e-9 

-3.03496e-9 

-3.03501e-9 

2 

20 Hz 

(1,1,0) 

ft 

-2.78908e-9 

-2.78875e-9 

-2.78917e-9 

-2.78921e-9 

3 

20 Hz 

(1.1,-40) 

A 

-5.14848e-ll 

-5.14711e-ll 

-5.14647e-ll 

-5.14996e-ll 

4 

20 Hz 

(1.1,-40) 

ft 

+4.97666e-10 

-|-4.97411e-10 

-|-4.97512e-10 

-|-4.97366e-10 

5 

40 Hz 

(1,1.0) 

A 

-1.07645e-9 

-1.07623e-9 

-1.07652e-9 

-1.07656e-9 

6 

40 Hz 

(1,1.0) 

ft 

-1.56155e-9 

-1.56153e-9 

-1.56178e-9 

-1.56180e-9 


We remark that the agreement between the two derivatives is remarkable, both for 
cases 1-4, where the wavefield is well-resolved, and for cases 5 and 6, where only 10 points 


objective functional, and not from the regularization part. Nevertheless, we have also successfully verified 
the accuracy of the regularization component of the gradients. 


20 









per wavelength are used for spatial discretization. 

The considered test problem is a heterogeneous half-space with a smoothly varying 
material profile along depth, given in (28), and mass density p = 2000 kg/m^, which, 
after truncation, is reduced to a cubic computational domain of length and width 24 m 
X 24 m, and 45 m depth. A 5 m-thick PML is placed at the truncation boundaries, as 
shown in Figure 5. The material properties at the interfaces T^ are extended into the PML. 
The interior and PML domains are discretized by quadratic hexahedral spectral elements 
(i.e., 27-noded bricks, and quadratic-quadratic pairs of approximation for displacement and 
stress components in the PML, and, also, quadratic approximation for material properties) 
of size 1 m, and At = 9 x 10“'^ s. Throughout, for the PML parameters, we choose 
Uo = 5, /3o = 400 and a quadratic profile for the attenuation functions, i.e., m = 2. 
See [13] for notation and other details. To probe the medium, we consider vertical stress 
loads with Gaussian pulse temporal signatures (see Table 1), applied on the surface of the 
domain over a region (—11 m < x,y < 11 m), whereas receivers that collect displacement 
response Um(x,t) are also located in the same region, at every grid point. To obtain 
synthetic data at the receiver locations, we use a model with identical characteristics and 
dimensions as detailed above, but, with a refined discretization; i.e., element size of 0.5 m, 
and At = 4.5 x 10“'^ s. The data was then mapped onto the coarser mesh discussed earlier. 
The total duration of the simulation is T = 0.5 s. 

4-2. Smoothly varying heterogeneous medium 

We consider a heterogeneous half-space with a smoothly varying material profile along 
depth, given by: 

/ (Izl — 22 51^\ 

A( 2 ) = p(z) = 80 + 0.45 \z\ -|- 35 exp (-^-) (MPa), (28) 

V 150 / 

and mass density p = 2000 kg/m^, which, after truncation, is reduced to a cubic compu¬ 
tational domain of length and width 40 m x 40 m, and 45 m depth. A 6.25 m-thick PML 
is placed at the truncation boundaries, as illustrated in Figure 6. The target profiles are 
shown in Fig, 7. The material properties at the interfaces F^ are extended into the PML. 
The interior and PML domains are discretized by quadratic hexahedral spectral elements 
(i.e., 27-noded bricks, and quadratic-quadratic pairs of approximation for displacement and 
stress components in the PML, and, also, quadratic approximation for material properties) 
of size 1.25 m, and At = 10“^ s. This leads to 3,578,136 state unknowns, and 616,850 
material parameters. To probe the medium, we consider vertical stress loads with Gaussian 
pulse temporal signatures (see Table 1), applied on the surface of the domain over a region 
(—17.5 m. < x,y < 17.5 m), whereas receivers that collect displacement response Um(x,t) 
are placed at every grid point, in the same region. 

Before attempting simultaneous inversion for the two Lame parameters, we perform 
single parameter inversion for a) p only, assuming A is a priori known and fixing it to the 
target profile; and b) A only, assuming distribution of p is known. 


21 




Figure 6: Problem configuration; material profile reconstruction of a smoothly varying medium. 



(a) 


(b) 


Figure 7: Smoothly varying medium: (a) target A and ji (MPa); and (b) profile at {x,y) = (0,0). 
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4-2.1. Single parameter inversion 

First, we assume A is a priori known, and fix it to the target profile. We start inverting 
for fi, with a homogeneous initial guess of 80 MPa, exploiting Tikhonov regularization for 
taming ill-posedness and solution multiplicity. We use the Gaussian pulse p 20 with maximal 
frequency content fmax = 20 Hz (see Table 1) for 50 iterations, and, then, switch to pao 
with fmax = 30 Hz. After 156 iterations, fi converges to the target profile, as shown in 
Figure 8(b). We compare the inverted cross-sectional prohles of with the target prohle at 
three different locations, as shown in Figure 8(c), 8(d), and 8(e). The agreement between 
the two prohles is excellent. Reduction of the misht functional with respect to inversion 
iterations is shown in Figure 10(a), which is almost 7 orders of magnitude. 

Next, we hx p, to the target prohle, and invert for A, starting with a homogeneous 
initial guess of 80 MPa. We use the Gaussian pulse p 2 o for 160 iterations, pso up to the 
300*^ iteration, and then switch to p 4 o. After 456 iterations, the optimizer converges to the 
prohle displayed in Figure 9(a). The agreement between the target prohle and the inverted 
prohle is remarkable. We compare the two prohles at three different cross-sections shown 
in Figure 9(c), 9(d), and 9(e): the agreement between the two prohles is excellent. The 
misht history is shown in Figure 10(b); the optimizer reduced the misht almost 6 orders of 
magnitude. 

We remark that the initial value of the misht in the hrst experiment is almost 2 orders 
of magnitude more than that of the second experiment. This indicates that the objective 
functional is not equally sensitive to both control parameters, as it has also been reported 
in [22]: the objective functional is more sensitive to p. 

4-2.2. Simultaneous inversion 

We start with a homogeneous initial guess of 80 MPa for both A and p and attempt 
simultaneous inversion. The target prohles are shown in Figure 7, and the inverted prohles 
are displayed in Figure 11(a) and 11(b). We also compare the cross-sectional values of 
the target and inverted prohles at three diherent locations, shown in Figure 12. Although 
the inverted p prohle agrees reasonably well with the target prohle, inversion for A is not 
satisfactory, and the inverted prohle departs from the target as depth increases. 

Due to the unsuccessful inversion of the A prohle in the case of simultaneous inversion, 
in the next experiment, we bias the search direction of A based on that of p, at the very 
early stages of inversion, according to the procedure detailed in Section 3.3.3. This leads to 
the successful reconstruction of the two prohles, as is shown in Figure 13(a) and 13(b). In 
Figure 14, we compare the cross-sectional values of the target and the inverted prohles. The 
agreement of the inverted p prohle with the target is remarkable. Moreover, the inverted A 
prohle agrees reasonably well with the target, with some discrepancies in depth. The misht 
history is shown in Figure 15(b), where the kink in the misht curve at the 50*^^ iteration 
corresponds to the termination point of the biasing scheme. 

We remark that in practical applications, one is more interested in the shear wave 
velocity (cs) and compression wave velocity (cp) prohles. Once the Lame parameters have 
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(a) A (a priori known) 


(b) /i (inverted) 





(c) p at {x,y) = (0,0) 


(d) p at (x,y) = (10,10) 


(e) p at {x,y) = (20, 20) 


Figure 8: Single-parameter inversion [p only) for a smoothly varying medium. 
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(a) A (inverted) 


(b) /i (a priori known) 





(c) A at {x,y) = (0,0) 


(d) A at {x, y) = (10,10) 


(e) A at {x,y) = (20,20) 


Figure 9: Single-parameter inversion (A only) for a smoothly varying medium. 
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misfit 




(a) reconstruction for /r only (b) reconstruction for A only 

Figure 10: Variation of the misfit functional with respect to inversion iterations (single parameter inversion). 



(a) A (inverted) (b) /i (inverted) 

Figure 11: Simultaneous inversion for A and /i using unbiased search directions. 
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(a) A at {x, y) = (0, 0) 


(b) A at {x,y) = (10,10) (c) A at {x,y) = (20,20) 



Figure 12: Cross-sectional profiles of A and y, using unbiased search directions. 
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(a) A (inverted) 


(b) ^ (inverted) 


Figure 13: Simultaneous inversion for A and fi using biased search directions. 


been determined, the wave velocities can be readily computed via: 







(29) 


In Figure 16, we compare the compression wave velocities at three different cross-sectional 
locations, where the agreement between the reconstructed Cp profile and the target is 
remarkable. The shear wave velocity does not depend on A, and, therefore, its quality is 
similar to that of the fj, profile. 


4-3. Layered medium 

We consider a 40 m x 40 m x 45 m layered medium, where a 6.25 m-thick PML is 
placed at its truncation boundaries. The properties of the medium are: 


'80 MPa, 


X{z) = ii{z) = < 


101.25 MPa, 
125 MPa, 


for — 12 m < z < 0 m, 
for — 27 m < 2 < —12 m, 
for — 50 m < 2 < —27 m. 


(30) 


and are shown in Figure 17, with mass density p = 2000 kg/m^. The material properties 
at the interfaces F^ are extended into the PML buffer. The interior and PML domains are 
discretized by quadratic hexahedral spectral elements (i.e., 27-noded bricks, and quadratic- 
quadratic pairs of approximation for displacement and stress components in the PML, and. 


28 









(a) A at {x, y) = (0, 0) 


(b) A at {x,y) = (10,10) (c) A at {x,y) = (20,20) 




-Target 

Inverted 

-Initial guess 




-Target 

Inverted 
- Initial guess 




-Target 

• Inverted 
- Initial guess 
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-15 
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-15 








■g 
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-30 
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-30 












(d) fi at {x, y) = (0, 0) (e) y, at {x, y) = (10,10) (f) y at {x, y) = (20, 20) 


Figure 14: Cross-sectional profiles of A and y using biased search directions. 
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misfit 




(a) unbiased search directions (b) biased search direction for A 

Figure 15: Variation of the misfit functional with respect to inversion iterations (simultaneous inversion). 





(a) Cp at {x, y) = (0, 0) 


(b) Cp at {x,y) = (10,10) (c) Cp at {x,y) = (20, 20) 


Figure 16: Cross-sectional profiles of Cp using biased search directions. 
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also, quadratic approximation for material properties) of size 1.25 m, and At = 10“^ s. For 
probing the medium, we use vertical stress loads with Gaussian pulse temporal signatures, 
applied on the surface of the domain over a region (—17.5m < x,y < 17.5m), whereas 
receivers that collect displacement response are also located in the same region, 

at every grid point. 



90 100 110 

A and /i (MPa) 


(a) 


(b) 


Figure 17: Layered medium: (a) target A and ^ (MPa); and (b) profile at {x,y) = (0,0). 

We start the inversion process with a homogeneous initial guess of 80 MPa for the 
Lame parameters, and attempt simultaneous inversion for both A and y, using the biasing 
scheme outlined in Section 3.3.3. We use the Total Variation regularization scheme, with 
e = 0.01, to capture the sharp interfaces of the target prohles. We use the Gaussian pulse 
P 20 , with fmax = 20 Hz, and hnal simulation time^^ T = 0.45 s, for 310 iterations. The 
resulting prohles are shown in Figure 18(a) and 18(b), where the layering of the medium 
is clearly visible in the inverted prohles. To improve the quality of the inverted prohles, 
we use them as an initial guess with the Gaussian pulse pso, and hnal simulation time of 
T = 0.4 s, for up to the 860**^ iteration, and, then, switch to p 4 o, with hnal simulation time 
of T = 0.4 s. After 1112 iterations, the optimizer converges to the prohles displayed in 


'^^The final simulation time is chosen such that they are long enough to probe the medium effectively, 
but not too long to increse the computational cost, unnecessarily. To this end, we run a forward simulation 
by using the target profile, and monitor when the total energy of the system (see [13]) die out. We then 
use this time duration for the inversion process. Reference [20] provides guidelines for choosing the source 
duration when inversion is performed by using field data. 


31 









Figure 18(c) and 18(d). There is excellent agreement between the inverted /i profile and 
the target profile. The inverted A profile is also in good agreement with the target profile: 
the two top layers have been reconstructed quite well, whereas the bottom layer is slightly 
“stiffer” in it’s middle zone. We compare the inverted profiles with the targets at three 
different cross-sections, shown in Figure 19. Due to the TV regularization, sharp interfaces 
have been captured quite successfully. In Figure 20, we compare the Cp profile with the 
target, at the same cross-sections; the agreement between the two profiles is impressive. 
Figure 21 shows the misfit history: the optimizer reduced the misfit almost 7 orders of 
magnitude. 


4.4- Layered medium with inclusion 

We consider a layered medium with an inclusion. The problem consists of a 40 m x 
40 m X 45 m layered medium with an ellipsoidal inclusion, where a 6.25 m-thick PML is 
placed at its truncation boundaries. The material profiles are given by: 


^80 MPa, 


X{z) = p( 2 ) = < 


101.25 MPa, 
125 MPa, 
156.8 MPa, 


for — 12 m < 2 < 0 m, 
for — 27 m < z < —12 m, 
for — 50 m < z < —27 m, 
for ellipsoidal inclusion, 


(31) 


and are shown in Figure 22, with constant mass density p = 2000 kg/m^. The ellipsoidal 
inclusion occupies the region ( ^~7-5 )2 _|_ _|_ (^z±n ^^2 ^ material properties at 

the interfaces F^ are extended into the PML buffer. The interior and PML domains are 
discretized by quadratic hexahedral spectral elements (i.e., 27-noded bricks, and quadratic- 
quadratic pairs of approximation for displacement and stress components in the PML, and, 
also, quadratic approximation for material properties) of size 1.25 m, and At = 10“^ s. To 
illuminate the domain, we use vertical stress loads with Gaussian pulse temporal signatures, 
applied on the surface of the medium over a region (—17.5m < x,y < 17.5m), whereas the 
receivers are also placed at every grid point in the same region. 

We use the Total Variation regularization scheme to alleviate ill-posedness and solution 
multiplicity, with e = 0.01. Similar to the previous examples, we use a source-frequency 
continuation scheme, starting with the Gaussian pulse p 2 o with maximal frequency content 
of 20 Hz for T = 0.45 s, and, when updates in the material profiles become practically 
insignificant, we switch to the next load in Table 1, which contains a broader range of 
frequencies, and, therefore, is able to image finer features. Figure 23(a) and 23(b) show the 
material profiles after 410 iterations, which adequately capture the layering of the domain 
as well as the ellipsoidal inclusion. To improve the quality of the reconstructed profiles, 
we use them as an initial guess with the Gaussian pulse pso, and final simulation time of 
T = 0.4 s, for up to the 730**^ iteration, and, then, switch to p 4 o, with final simulation time 
of r = 0.4 s. Figure 23(c) and 23(d) show the inverted profiles after 1160 iterations. The 
sharp interfaces between the three layers and around the ellipsoidal inclusion are very well 
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(c) A {fmax = 40 Hz) (d) fj, {fmax = 40 Hz) 


Figure 18: Simultaneous inversion for A and /r (layered medium). 
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A (MPa) A (MPa) A (MPa) 


(a) A at {x, y) = (0, 0) 


(b) A at {x,y) = (10,10) (c) A at {x,y) = (20,20) 



fi (MPa) y. (MPa) y (MPa) 


(d) fi at (x, y) = (0, 0) 


(e) y at {x,y) = (10,10) (f) y at {x,y) = (20, 20) 


Figure 19: Cross-sectional profiles of A and y (layered medium). 
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(a) Cp at {x,y) = (0,0) (b) Cp at {x,y) = (10,10) (c) Cp at {x,y) = (20, 20) 


Figure 20: Cross-sectional profiles of Cp (layered medium). 



Figure 21: Variation of the misfit functional with respect to inversion iterations (layered medium). 
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(a) 


(b) 


Figure 22: Layered medium with inclusion: (a) target A and /r; and (b) profile at {x,y) = (7.5,0). 

captured for the p profile. The A profile agrees reasonably well with the target, showing 
some “stiff’ features at the center of the bottom layer, similar to the previous example. 

Figures 24 and 25 compare the inverted profiles with the target profiles at three different 
cross-sectional lines of the domain, indicating successful imaging of both the layering and 
the inclusion. Variation of the misfit functional with respect to the inversion iterations 
is shown in Figure 25, where, again, a kink at the 50*^^ iteration of the misfit curve, 
corresponds to the termination point of the biasing scheme. 

Encouraged by the successful performance of the proposed inversion algorithm with 
noise-free data, next, we consider adding different levels of Gaussian noise to the measured 
synthetic response at the receiver locations, and investigate its effect on the inversion. 
Figures 27(a)-27(d) show the measured displacement response of the system at {x,y,z) = 
(3.125,-13.75,0) m, subjected to the p 2 o pulse, contaminated with 1%, 5%, 10%, and 
20% Gaussian noise, respectively. Using the source-frequency continuation scheme, the 
optimizer converges after 811 and 751 iterations, respectively, for cases corresponding to 
the 1% and 5% Gaussian noise levels. The inverted profiles are shown in Figure 28. The 
reconstruction is successful, with minor discrepancies on the top surface. Next, we increase 
the noise level to 10% and 20%, and attempt inversion; after 770 and 674 iterations, 
respectively, we converge to the profiles shown in Figure 29. The quality of the inverted 
prohles decreases as the noise level increases. However, similarly to the previous case, 
except for a thin layer on the top surface, inversion is successful. In Figure 30, we compare 
cross-sectional profiles of A and p with the target, at different noise levels, at {x, y) = 
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(a) A {fmax = 10 Hz) 


(b) fj. {fmax = 10 Hz) 



(c) A {fmax = 40 Hz) 


(d) fj, (fmax = 40 Hz) 


Figure 23: Simultaneous inversion for A and ^ (layered medium with inclusion). 
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(a) A at {x,y) = (7.5,0) (b) A at {x,y} = (7.5,10) (c) A at {x,y) = (20,20) 



(d) fi at {x,y) = (7.5,0) (e) y. at {x,y) = (7.5,10) (f) y at {x,y) = (20,20) 

Figure 24: Cross-sectional profiles of A and y (layered medium with inclusion). 
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c, (a/,) c, (m/.) c, (n/.) 


(a) Cp at (x,y) = (7.5,0) (b) Cp at (x,y) = (7.5,10) (c) Cp at (x,y) = (20, 20) 


Figure 25: Cross-sectional profiles of Cp (layered medium with inclusion). 



Figure 26: Variation of the misfit functional with respect to inversion iterations (layered medium with 
inclusion). 
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(7.5,0) m, which passes through the center of the ellipsoidal inclusion. Sharp interfaces 
are captured remarkably well for the fi prohle, even at the presence of 20% noise. The 
inversion for A is also satisfactory. 



(a) with 1% Gaussian noise 



(c) with 10% Gaussian noise 



(b) with 5% Gaussian noise 



(d) with 20% Gaussian noise 


Figure 27: Measured displacement response of the layered medium with inclusion, at {x,y,z) = 
(3.125, —13.75, 0) m, due to the Gaussian pulse p 2 o, contaminated with Gaussian noise. 

4-5. Layered medium with three inclusions 

In the last example, we consider a layered medium, with three inclusions, to study the 
performance of our inversion scheme for a more complex material profile. The problem 
consists of an 80 m X 80 m x 45 m medium, where a 6.25 m-thick PML is placed at its 
truncation boundaries. The material prohles are given by: 
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60 170 

(c) A (5% Gaussian noise) (d) ji (5% Gaussian noise) 


Figure 28: Simultaneous inversion for A and /i using measured data contaminated with 1% and 5% Gaussian 
noise (layered medium with inclusion). 
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(c) A (20% Gaussian noise) (d) fi (20% Gaussian noise) 


Figure 29: Simultaneous inversion for A and fi using measured data contaminated with 10% and 20% 
Gaussian noise (layered medium with inclusion). 
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(a) A at {x,y) = (7.5,0) 


(b) 11 at {x,y) = (7.5,0) 


Figure 30; Cross-sectional profiles of A and y at different noise levels (layered medium with inclusion). 


for — 15 m < 2 < 0 m, 
for — 30 m < 2 ; < —15 m, 
for — 50 m < 2 ; < —30 m, 

for spheroidal: < h 

for ellipsoidal: + (^)2 < 1, 

for sphere: {x — 20)^ + (j/ + 20)^ + (2 + 35)^ < 6.25, 

and are shown in Figure 31, with constant mass density p = 2000 kg/m^. Figures 32(a) 
and 32(b) depict the target profiles on a cross-section through the domain situated at 
8.75 m and 35 m from the top surface, going through the ellipsoid’s and sphere’s midplane, 
respectively. In terms of the smallest wavelength^^ the prescribed geometry comprises a 
domain of 16 x 16 x 9 wavelengths long, wide, and deep, a spherical inclusion with a diameter 
2.5 wavelengths, an ellipsoidal inclusion of 6 x 3 x 2 wavelengths, and a spheroidal inclusion 
of 1.5 X 8 X 1.5 wavelengths. The material properties at the interfaces F^ are extended into 
the PML buffer. The interior and PML domains are discretized by quadratic hexahedral 


A( 2 ) = p.{z) = < 


80 MPa, 
101.25 MPa, 
125 MPa, 
156.8 MPa, 
156.8 MPa, 
80 MPa, 


^^The smallest wavelength is equal to the smallest velocity in the formation 200 m/s, divided by the 
largest probing frequency 40 Hz, i.e., 5 m. 
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spectral elements (i.e., 27-noded bricks, and quadratic-quadratic pairs of approximation 
for displacement and stress components in the PML, and, also, quadratic approximation 
for material properties) of size 1.25 m, and At = 10“^ s. This leads to 9,404,184 state 
unknowns, and 2,429, 586 material parameters. To illuminate the domain, we use vertical 
stress loads with Gaussian pulse temporal signatures, applied on the surface of the medium 
over a region (—37.5 m < x,y < 37.5 m), whereas receivers are placed at every grid point, 
within the same region as the load. 

To narrow the feasibility space and alleviate difficulties with solution multiplicity, we 
use the Total Variation regularization, with e = 0.01, combined with the regularization 
factor continuation scheme outlined in Section 3.3.1, the source-frequency continuation 
scheme in Section 3.3.2, and the biasing scheme for A search directions in Section 3.3.3. 
Specifically, we use the regularization parameter g = 0.5 when illuminating the medium 
with pulse P 20 for 60 iterations. Next, we use g = 0.4 with pulse pso up to the 290*'^ 
iteration. Finally, we use g = 0.3 with pulse p 4 o and stop at the 741®* iteration. In all the 
three cases, the total simulation time is T = 0.7 s. 

Figure 33 shows the inverted prohle along a cross-section that cuts through the domain 
from {x,y) = (-20,-46.5) to (—20,20) to (46.5,20). The cross section passes through 
the larger semi-principal axes of both ellipsoids, and shows very good reconstruction of 
the y profile, and satisfactory inversion of the A prohle. The layering is sharp, and the 
ellipsoids are captured well. In Figure 34, a cross section of the inverted prohles from 
{x, y) = (20,46.5) to (20, —20) to (—46.5, —20) is displayed, where it passes through the 
smaller semi-principal axes of the ellipsoids and the center of the sphere. The ellipsoids 
are well captured; however, the sphere, which consists of “soft” materials, can hardly be 
noticed, especially, in the A prohle. Figure 35(a) and 35(b) show the inverted prohles on a 
cross-section through the domain, situated at 8.75 m from the surface, going through the 
top ellipsoid’s midplane, and show satisfactory reconstruction of the ellipsoid. To see the 
reconstruction of the sphere in more detail, we consider a cross-section, which goes through 
the sphere’s midplane, situated at 35 m from the top surface; this is shown in Figure 35(c) 
and 35(d). The sphere’s footprint is visible in the A prohle, whereas it is more conspicuous 
in the y prohle. 

We also compare cross sections of the inverted prohles with the target along three 
different lines, which pass through the ellipsoids and the sphere. These are shown in 
Figure 36. Overall, the inverted prohles are satisfactory. 

5. Conclusions 

We discussed a full-waveform-based inversion methodology for imaging the elastic prop¬ 
erties of a soil medium in three-dimensional, arbitrarily heterogeneous, semi-inhnite do¬ 
mains. The problem typically arises in geotechnical site characterization and geophysical 
explorations, where high-hdelity imaging of the two Lame parameters (or an equivalent 
pair) is of interest. Elastic waves are used as probing agents to interrogate the soil medium. 
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(b) 



A and (MPa) A and ^ (MPa) A and fi (MPa) 


(c) 


(d) 


Figure 31: Layered medium with three inclusions: target A and fi (a) along a cross-section that cuts through 
the domain from {x,y) = (-20,-46.5) to (-20,20) to (46.5,20); (b) along a cross-section that cuts through 
the medium from {x,y) = (20,46.5) to (20,-20) to (^6.5,-20); (c) profile at {x,y) — (-20,-20); (d) profile 
at {x,y) = (20,20); and (e) profile at {x,y) = (20,-20). 
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Figure 32: Layered medium with three inclusions: target A and /r on (a) the z = —8.75 m cross-section; 
and (b) the z = —35 m cross-section. 




(a) A {fmax = 40 Hz) 


(b) ^ [fmax = 40 Hz) 


Figure 33: Simultaneous inversion for A and fi'. cross-section cuts through the domain from {x, y) = 
(-20,-46.5) to (—20,20) to (46.5,20) (layered medium with three inclusions). 
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(a) A {fmax = 40 Hz) 


(b) {fmax = 40 Hz) 


Figure 34: Simultaneous inversion for A and /r: cross-section cuts through the domain from (k, y) = (20, 46.5) 
to (20, —20) to (—46.5, —20) (layered medium with three inclusions). 


and the response of the medium to these waves are collected at receivers located on the 
ground surface. The inversion process relies on minimizing a misht between the collected 
response at receiver locations, and a computed response based on a trial distribution of 
the Lame parameters. We used the apparatus of PDE-constrained optimization to im¬ 
pose the forward wave propagation equations to the minimization problem, directly in the 
time-domain. Moreover, PMLs were used to limit the extent of the computational domain. 

To alleviate the ill-posedness, associated with inverse problems, we used regularization 
schemes, along with a regularization factor continuation scheme, which tunes the regular¬ 
ization factor adaptively at each inversion iteration. We discussed additional strategies to 
robustify the inversion algorithm: specifically, we used (a) a source-frequency continuation 
scheme such that the inversion process evolves by using low-frequency sources, and, grad¬ 
ually, we use sources with higher frequencies; and (b) a biasing scheme for the A-prohle, 
such that, at early iterations of inversion, the search direction for A is biased based on that 
of /r. The latter strategy, in particular, improves the reconstruction of the material profiles 
when simultaneous inversion of the two Lame parameters is exercised. To the best of our 
knowledge, this is the hrst attempt that the two Lame parameters have been successfully 
reconstructed in three-dimensional PML-truncated domains. 

In order to resolve the forward wave propagation problem, we used a recently developed 
hybrid hnite element approach, where a displacement-stress formulation for the PML is 
coupled to a standard displacement-only formulation for the interior domain, resulting in 
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(a) A {fmax = 40 Hz) 


(b) {fmax = 40 Hz) 




(c) A {fmax = 40 Hz) 


(d) II {fmax = 40 Hz) 


Figure 35: Layered medium with three inclusions: (a) inverted A profile on the z = —8.75 m cross-section; 
(b) inverted /r profile on the 2 : = —8.75 m cross-section; (c) inverted A profile on the 2 = —35 m cross-section; 
and (d) inverted /r profile on the 2 = —35 m cross-section. 
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(a) A at (x, y) = (—20, —20) 



f (MP.) 


(d) y at [x, y) = (-20, -20) 



A (MPa) 


(b) A at {x, y) = (20,20) 



y. (MPa) 

(e) y at {x,y) = (20,20) 



A (MPa) 


(c) A at {x,y) = (20,-20) 



p (MP.) 

(f) y at {x,y) = (20, -20) 


Figure 36: Cross-sectional profiles of A and y (layered medium with three inclusions). 
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a scheme with optimal computational cost. Time-integration is accomplished by using 
an explicit Runge-Kutta scheme, which is well-suited for large-scale problems on parallel 
computers. 

By comparing directional finite differences of the discrete objective functional, and 
directional derivatives obtained via the control problems, we verihed the accuracy and 
correctness of the material gradients. We reported numerical results demonstrating suc¬ 
cessful reconstruction of both Lame parameters for smooth and sharp profiles. Overall, the 
framework discussed in this paper seems practical, and promising. 
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Appendix A. Gradient of a fnnctional 

The gradient of a functional J- : H ^ where 77 is a Hilbert space, is defined as the 
Riesz-representation of the derivative .T'(q)(q), such that: 

( 6^(q),q )h =-^'(q)(q) Vqe77, (A.i) 

where G denotes the gradient, and we use the following notation for the Gateaux derivative 
of at q in a direction q: 


(A^2) 

/i->0 n 

With this definition, it is not possible to talk about the gradient, without specifying the 
inner-product utilized to represent the derivative [5]. 

Appendix B. The adjoint problem time-integration scheme 

We outline the explicit T^-order Runge-Kutta method (RK-4) for the reverse time- 
integration of the adjoint problem. Upon using spectral elements, with Legendre-Gauss- 
Lobatto (LGL) quadrature rule, the mass-like matrix M becomes diagonal; therefore, its 
inverse can be readily computed. We use the following notation: 
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C = C M-\ 
G = G 


-1 


Using (B.l), (15) becomes: 

dt 


Vo' 

yi 

= 

ya. 



0 

0 


I 

0 


0 

I 

CT 


LG^ -K^ 

The scheme entails computing the following vectors 


K = K M 

f = M-i 



'yo' 


o' 


yi 

+ 

0 


y2. 


f_ 


(B.la) 

(B.lb) 


(B.2) 


kio = y?, 
kii = y2, 

ki2 = Cy^ - Ky^ + Gy^ + t, 
k2o = yi - ^ kii, 
k21 = y2 - ^ ki2, 

k 22 = C(y^ - ^ ki 2 ) - K{yi - ^ kii) + G(y^ - ^ kio) + 
kso = yi - ^ k 2 i, 

ksi = y2 - ^ k22, 

k32 = C(y^ - ^ k22) - K{yi - Y ^2i) + G(y^ - ^ k20) + 
k40 = yi - At k3i, 

k41 = y2 “ At k32, 

k42 = C(y^ - At k32) - K(y5^ - At k3i) + G(y^ - At k3o) + P'^ 

Finally, the solution at time step (n — 1) can be updated via: 


'yo' 

n—1 

'yo 

n 

At 

kio + 2 k2o + 2 k3o + kio 

yi 

= 

yi 


kii + 2 k2i + 2 koi + kii 

y2. 


y2. 

D 

ki2 + 2 k22 + 2 k32 + k42 


(B.4) 
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Appendix C. Discretization of the control problems 

In Section 3.1.3, we discussed the A- and /i-control problems. In this part, we consider 
their spatial discretization. As discussed in [13], we use the basis $ for the spatial dis¬ 
cretization of w(x, t) and n(x, t), and x is the basis for discretizing A(x) and /i(x). For 
instance, if we approximate A(x) with A/i(x), then A/i(x) = x^A, where A comprises the 
vector of nodal values for A. In the following, subscripts in the basis indicate derivatives, 
and u/j = (u^,ny ,nj)^ and Wh = (w!^, wj, is the vector of discrete values of the 

state and adjoint variables, respectively. Accordingly: 


M = XX^ dfl. 

For Tikhonov regularization: 

Sreg = jiXxXx + XyXy + XzXl)>^ dO, 

Sreg = jiXxXx + XyXy + XzXl)tJ- dfl. 


For Total Variation regularization: 

iXxXl + XyXy + XzXz)^^ 


= 

Sreg 




= 

&reg 




[^'^{XxXx + XyXy + XzXT)>^ + e 
iXxXl + XyXy + XzXz)tJ' 
[fJ-'^iXxXx + XyXy + XzX^)l^ + e 


dfl. 


dn. 


(C.la) 


(C.lb) 


(C.lc) 


(C.ld) 


(C.le) 


Moreover, 

gmis = -y j X {^l^x + ^yy^y + ^'^^z){^lu:, + ^'^Uy + ^'^U;,)dn dt, (C.lf) 

Smis = “ y y ^ + ^yWy ^'^Uy + W, 

+ {^yV^X + ^lWy){^lUy + ^yU^) T ^ X T U,, ) 

Uj,)^ dfi dt. (C.lg) 

In (C.la), upon using spectral elements with LGL quadrature rule, M becomes diagonal; 
thus, its inverse can be computed easily. 
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